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Abstract — We present a new algorithm for solving a poly- 
nomial program P based on the recent "joint + marginal" 
approach of the first author for parametric polynomial opti- 
mization. The idea is to first consider the variable Xi as a 
parameter and solve the associated (n — 1) -variable (x2, . . . , x n ) 
problem P(ki) where the parameter xi is fixed and takes values 
in some interval Yi C R, with some probability tpi uniformly 
distributed on Yi. Then one considers the hierarchy of what 
we call "joint+marginal" semidefinite relaxations, whose duals 
provide a sequence of univariate polynomial approximations 
xi n> pk(xi) that converges to the optimal value function 
xi n> J{x\) of problem P(xi), as k increases. Then with k fixed 
a priori, one computes x\ € Yi which minimizes the univariate 
polynomial Pk{xi) on the interval Yi, a convex optimization 
problem that can be solved via a single semidefinite program. 
The quality of the approximation depends on how large k can 
be chosen (in general for significant size problems k = 1 is the 
only choice). One iterates the procedure with now an (n — 2)- 
variable problem P(a;2) with parameter X2 in some new interval 
Y2 C R, etc. so as to finally obtain a vector x £ R™. Preliminary 
numerical results are provided. 

I. Introduction 



Consider the general polynomial program 

P : /* := min{/(x) : x G K } 



(1) 



where / is a polynomial, K C K™ is a basic semi-algebraic 
set, and /* is the global minimum of P (as opposed to a local 
minimum). One way to approximate the global optimum 
/* of P is to solve a hierarchy of either LP-relaxations or 
semidefinite relaxations as proposed in e.g. Lasserre [4], [5]. 
Despite practice with the semidefinite relaxations seems to 
reveals that convergence is fast, the matrix size in the z-th 
semidefinite relaxation of the hierarchy grows up as fast as 
0(n l ). Hence, for large size (and sometimes even medium 
size) problems, only a few relaxations of the hierarchy can 
be implemented (the first, second or third relaxation). In that 
case, one only obtains a lower bound on /*, and no feasible 
solution in general. So an important issue is: 

How can we use the result of the i-th semidefinite relax- 
ation to find an approximate feasible solution of the original 
problem? 

For some well-known special cases of 0/1 optimization 
like e.g. the celebrated MAXCUT problem, one may gener- 
ate a feasible solution with guaranteed performance, from a 
randomized rounding procedure that uses an optimal solution 
of the first semidefinite relaxation (i.e. with i — 1); see 
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Goemans and Williamson [2]. But in general there is no 
such procedure. 

Our contribution is to provide two relatively simple al- 
gorithms for polynomial programs which builds up upon 
the so-caUed "joint+marginal" approach (in short (J+M)) 
developed in [6] for parametric polynomial optimization. 
The (J+M)-approach for variables x G R" and parameters 
y in a simple set Y, consists of the standard hierarchy of 
semidefinite relaxations in [4] where one treats the param- 
eters y also as variables. But now the moment-approach 
implemented in the semidefinite relaxations, considers a joint 
probability distribution on the pair (x, y), with the additional 
constraint that the marginal distribution on Y is fixed (e.g. 
the uniform probability distribution on Y); whence the name 
"joint+marginal " . 

For every k = l,...,n, let the compact interval Y& := 
[x fc ,Xfe] C t be contained in the projection of K into the 
Xfc-coordinate axis. In the context of the (non-parametric) 
polynomial optimization (fTJ, the above (J+M)-approach can 
be used as follows in what we call the (J+M)-algorithm: 

• (a) Treat x\ as a parameter in the compact interval Yi = 
[•Ei> w hh associated probability distribution ipi uniformly 
distributed on Yi. 

• (b) with i € N fixed, solve the i-th semidefinite 
relaxation of the (J+M)-hierarchy [6] applied to problem 
P(:ei) with n — 1 variables (#2, ■ • ■ > x n ) and parameter x\, 
which is problem P with the additional constraint that the 
variable xi G Yi is fixed. The dual provides a univariate 
polynomial x% H> J}{xi) which, if i would increase, would 
converge to J 1 (xi) in the Li(< / 9i)-norm. (The map v M> 
J 1 (w) denotes the optimal value function of P(w), i.e. the 
optimal value of P given that the variable x\ is fixed at the 
value v.) Next, compute x% G Yi, a global minimizer of 
the univariate polynomial J\ on Yi (e.g. this can be done 
by solving a single semidefinite program). Ideally, when i is 
large enough, x% should be close to the first coordinate x^of 
a global minimizer x* = x*, . . . , a;*) of P. 

• (c) go back to step (b) with now G Y2 C M instead 
of X\, and with ip 2 being the probability measure uniformly 
disttibuted on Y 2 . With the same method, compute a global 
minimizer i 2 £ Y 2 , of the univariate polynomial x 2 *-> 
Jj 2 (x 2 ) on the interval Y 2 . Again, if i would increase, Jf 
would converge in the Li(<^ 2 )-norm to the optimal value 
function v H» J 2 (v) of P(x 2 ) (i.e. the optimal value of P 
given that the variable x 2 is fixed at the value v.) Iterate until 
one has obtained x n G Y„ c K. 

One ends up wih a point x G rjfc=i "^k and in general 
x ^ K. One may then use x as initial guess of a local 
optimization procedure to find a local minimum x G K. 



The rational behind the (J+M)-algorithm is that if i is large 
enough and P has a unique global minimizer x* G K, then 
x as well as x should be close to x*. 

The computational complexity before the local optimiza- 
tion procedure is less than solving n times the i-th semidef- 
inite relaxation in the (J+M)-hierarchy (which is itself of 
same order as the i-th semidefinite relaxation in the hierarchy 
defined in [4]), i.e., a polynomial in the input size of P. 

When the feasible set K is convex, one may define the 
following variant to obtain a feasible point x G K. Again, 
let Yi be the projection of Ki into the x i -coordinate axis. 
Once xx £ Yi is obtained in step (b), consider the new opti- 
mization problem P(ii) in the n—1 variables (22, • • • , x n ), 
obtained from P by fixing the variable x\ G Yi at the 
value x\. Its feasible set is the convex set Ki := K n {x : 
%l = xi}. Let Y2 be the projection of Ki into the in- 
coordinate axis. Then go back to step (b) with now X2 G Y 2 
as parameter and (x^, . . . , x n ) as variables, to obtain a point 
$2 S Y2, etc. until a point x G IIfc=i ^fe * s obtained. Notice 
that now x € K because K is convex. Then proceed as 
before with x being the initial guess of a local minimization 
algorithm to obtain a local minimizer x € K of P. 

II. The "joint+marginal approach to parametric 

OPTIMIZATION 

Most of the material of this section is taken from [6], 
Let K[x, y] denote the ring of polynomials in the variables 
x = (xx,...,x n ), and the variables y = (yi,...,y p ), 
whereas R[x, y]d denotes its subspace of polynomials of 
degree at most d. Let S[x, y] C R[x, y] denote the subset of 
polynomials that are sums of squares (in short s.o.s.). For a 
real symmetric matrix A the notation A y stands for A 
is positive semidefinite. 

The parametric optimization problem 

Let Y C W be a compact set, called the parameter set, 
and let /, hj e M[x], j = 1, . . . , m. Let K c R™ x W be 
the basic closed semi-algberaic set: 

K:={(x,y) : y G Y ; fy(x,y) > 0, j = l,...,m} (2) 

and for each y G Y, let 

K y := {xGR" : (x, y) G K}. (3) 

For each y G Y, fixed, consider the optimization problem: 

J(y) := inf{/(x,y) : (x,y) G K}. (4) 

X 

The interpretation is as follows: Y is a set of parameters 
and for each instance y G Y of the parameter, one wishes 
to compute an optimal decision vector x*(y) that solves 
problem (@). Let ip be a Borel probability measure on Y, with 
a positive density with respect to the Lebesgue measure on 
M. p (or with respect to the counting measure if Y is discrete). 
For instance 

ip(B) := ( [ dy) f dy, MB e B(W), 
\Jy J JYnB 



is uniformly distributed on Y. Sometimes, e.g. in the context 
of optimization with data uncertainty, ip is already specified. 
The idea is to use <p (or more precisely, its moments) to get 
information on the distribution of optimal solutions x*(y) 
of P y , viewed as random vectors. In this section we assume 
that for every y G Y, the set K y in © is nonempty. 

A. A related infinite-dimensional linear program 

Let M(K) be the set of finite Borel probability measures 
on K, and consider the following infinite-dimensional linear 
program P: 

P ■= inf \ f dp, ; np, = <p \ , (5) 

C6M(K) IJ K J 

where nfj, denotes the marginal of /1 on R p , that is, irp, is a 
probability measure on R p defined by nfi(B) := /i(R" x B) 
for all B G B(W). Notice that /i(K) = 1 for any feasible 
solution /i of P. Indeed, as ip is a probability measure and 
7r/i = ip one has 1 = ip(Y) = fi(R n x W) = p,(K). 

The dual of P is the the following infinite-dimensional 
linear program: 

P* ■= sup / p(y)d(p(y) 

p£K[y] JY (O) 

/(x)-p(y) > V(x,y) GK. 

Recall that a sequence of measurable functions (g n ) on 
a measure space (Y , B(Y), ip) converges to g, ip-almost 
uniformly, if and only if for every e > 0, there is a set 
A G B(Y) such that (p(A) < e and g n — > g, uniformly on 
Y\A. 

Theorem 1 ([6]): Let both Y C W and K in © be 
compact and assume that for every y G Y, the set K y C R" 
in (01 is nonempty. Let P be the optimization problem © 
and let X* := {x G R" : /(x,y) = J(y)}, y G Y. Then: 

(a) p = J J(y) dip(y) and P has an optimal solution. 

(b) Assume that for (^-almost y G Y, the set of minimizers 
of X* is the singleton {x*(y)} for some x*(y) G K y . Then 
there is a measurable mapping g : Y — > K y such that 

g(y) = x*(y) for every y G Y 

p = [/(s(y)j)My), (7) 

and for every a G N", and (3 eW: 

f xYf(xj) = / y%(y) tt ^(y). (8) 

(c) There is no duality gap between © and ©, i.e. p = p* , 
and if (pj)ieN C R[y] is a maximizing sequence of © then: 

J \J{y) -Pi(y)\d<p(y) asi^oo. (9) 
Moreover, define the functions (pi) as follows: po := po, and 

y-^Pi(y) := max[pi_i(y),pj(y)], i = l,2,... 
Then pi — > J(-), 93-almost uniformly. 



An optimal solution p* of P encodes all information on 
the optimal solutions x*(y) of P y . For instance, let B be a 
given Borel set of R™ . Then from Theorem Q] 

Prob(x*(y) G B) = p*(B x R p ) = <p{g- l {B)), 

with g as in Theorem |TJb). 

Moreover from Theorem \V[c), any optimal or nearly 
optimal solution of P* provides us with some polynomial 
lower approximation of the optimal value function y h-> J(y) 
that converges to «/(■) in the L\(<p) norm. Moreover, one 
may also obtain a piecewise polynomial approximation that 
converges to J(-), ^-almost uniformly. 

In [6] the first author has defined a (J+M)-hierarchy of 
semidefinite relaxations (Q,) to approximate as closely as 
desired the optimal value p. In particular, the dual of each 
semidefinite relaxation provides a polynomial g< G E[y] 
bounded above by J(y), andy v+ q. t (y) := max^ = i,.„j <#(y) 
converges ip-almost uniformly to the optimal value function 
J, as i — > oo. This last property is the rationale behind the 
heuristic developed below. 

III. A "JOINT+MARGINAL" APPROACH 

Let := {a G N n : |a| < i} with |a| = Yn<*i- With 
a sequence z = (z a ) indexed in the canonical basis (x a ) of 
R[x], let L z : R[x] — > K be the linear mapping: 

/(=£/ a (x)) ^ L M (/) := J2f-z a , f G K[x]. 

a a 

Moment matrix: The moment matrix M,(z) associated 
with a sequence z = (z Q ), a G NJj, has its rows and columns 
indexed in the canonical basis (x Q ), and with entries. 

Mi(*)(a,0) = L z (x a+ ?) = z a+0 , Ma, (3 G N?. 

Localizing matrix: Let q be the polynomial x i— > q(x) :— 
^2 U q u X- u . The localizing matrix Mj(g z) associated with q G 
R[x] and a sequence z = has its rows and columns 

indexed in the canonical basis (x Q ), and with entries. 

Mi(9z)(a,/3) = L.(<?(x)x*^) 

= <?«2:q+/J+« 5 Va,/3eN". 

A sequence z = (z Q ) C 1 is said to have a representing 
finite Borel measure supported on K if there exists a finite 
Borel measure p such that 

z a = [ * a dp, Va G N™. 

A. A "joint+marginal" approach 

With C M[x], let K C R" be the basic 

compact semi-algebraic set 

K := {x e R" : &(x) > 0, J = 1, . . . , m}, (10) 

and consider the polynomial optimization problem ((T). 

Let Yfc C R be some interval [x k ,Xk], assumed to be 
contained in the orthogonal projection of K into the x k - 
ccordinate axis. 



For instance when the gj's are affine (so that K is a 
convex polytope), x k (resp. Wk) solves the linear program 
min(resp max) {xk : x G K}. Similarly, when K is convex 
and defined by concave polynomials, one may obtain x k and 
Xk, up to (arbitrary) fixed precision. In many cases, (upper 
and lower) bound constraints on the variables are already 
part of the problem definition. 

Let ifik the probability measure uniformly distributed on 
Yfc, hence with moments {fit) given by: 

(3 e = / x k d Vk {x) = ( *- k ; (ID 
Jx, {k + l){x k -x k ) 

for every £ = 0,1,.... Define the following parametric 
polynomial program in n — 1 variables: 

J k {y) = min{/(x) : x G K; x k = y}, (12) 

X 

or, equivalently J k (y) — min{/(x) : x G K^}, where for 
every y G Y: 

K y := {x G K; x k = y}. (13) 
Observe that by definition, /* = min{J fc (x) : x G Y^}, 

X 

and K. y ^ whenever y G Y k , where Y^ is the orthogonal 
projection of K into the x k -coordinate axis. 

Semidefinite relaxations 

To compute (or at least approximate) the optimal value 
p of problem P in (0) associated with the parametric 
optimization problem (TT2b . we now provide a hierarchy of 
semidefinite relaxations in the spirit of those defined in [4]. 

Let Vj := [(deggj)/2~|, j = 1, . . . ,m, and for i > maxj Vj, 
consider the semidefinite program: 

p tk = inf L z {f) (14) 

Z 

s.t. Mi(z) h 0, Mi-*, {gj z) t 0, j = 1, . . . ,m 
L z {xi)=/3 £ , £ = 0,l,...2i, 

where (/3^) is defined in (fTTT i. We call dT~4l> the parametric 
semidefinite relaxation of P with parameter y — x k . Observe 
that without the "moment" constraints L z {x k ) — fy, £ = 
1, . . . 2i, the semidefinite program ( TT~4T > is a relaxation of P 
and if K is compact, its corresponding optimal value /* 
converges to /* as k oo; see Lasserre [4]. 
Letting 90 = 0, the dual of (TBI) reads: 

Pik = SU P Xe 

2i tci 

s-t /(x)-^Af4=a + ^a jft (15) 
fcO 3=1 
(jj G S[x], < j < m; 
degajgj <2i, < j < m. 

Equivalently, recall that R^^i is the space of univariate 
polynomials of degree at most 2i, and observe that in dl~5b . 
the criterion reads 

^\ i p i = I Pi (y)d(p k (y), 
£=0 Yfc 



where pi 6 R[xfc]2i is the univariate polynomial x% h4 
Pi{xk) '■= YlgLo^e x k- Then equivalently, the above dual 
may be rewritten as: 

p* k = sup / pidifk 

m 

S.t. / - ft = <T + 53 ^ ^ 16 ) 

3=1 

K S M[a; fc ]2i; o-j G £[x], < j < to; 
degOjgj <2i, < j < m. 

Assumption 1: The family of polynomials (g^) C R[x] is 
such that for some M > 0, 

m 

x^M-\\xf = a +J2 (X 3 9j, 

3=1 

for some M and some s.o.s. polynomials (<7j) C E[x]. 

Theorem 2: Let K be as ( fTQb and Assumption [TJ hold. 
Let the interval Yfe C 1 be the orthognal projection of K 
into the Xfe-coordinate axis, and let <p k be the probability 
measure, uniformly distributed on Yfe. Assume that J£ y in 
(TT~3T > is not empty, let y h-> J fc (y) be as in ( TT~2T > and consider 
the semidefinite relaxations (TT4T> - (fT~6b . Then as i — > oo: 

(a) p ife t / J k d<pk and p* fe t / J k dip k 

(b) Let (pi, ((Tj)) be a nearly optimal solution of ( fToT l, e.g. 
such that J Yk Pidip k > p* k - Then < J k (y) for 
all y G Yfe, and 

/ l^ fc (y) -Pi(y)\difk(y) -> 0, as i -> oo. (17) 

Moreover, if one defines po := Po> an d 

y-^Pt{y) ■= ioax[pi-i{y),Pi(y)], i = i,2,..., 

then t J k (y), for <y9fc-almost all y € Yfe, and so pi — > 
J k , iy9fe-almost uniformly on Y^. 

Theorem |2] is a direct consequence of [6, Corollary 2.6]. 

B. A "joint+marginal" algorithm for the general case 

Theorem [2] provides a rationale for the following (J+M)- 
algorithm in the general case. In what follows we use the 
primal and dual semidefinite relaxations (TT~4-b - (fT3T> with index 
i fixed. 

ALGO 1: (J+M)-algorithm: non convex K, relaxation i 

Set jfe = 1; 

Step k: Input: K, /, and the orthogonal projection Y^ = 
[•Efci^fe] °f K into the Xk -coordinate axis, with associated 
probability measure ipk, uniformly distributed on Yfe. 
Ouput: x k e Y fe . 

Solve the semidefinite program ( fToT l and from an optimal 
(or nearly optimal) solution (j>i,{o~j)) of (TToT l. get a global 
minimizer x k of the univariate polynomial pi on Yfe. 
If k = n stop and output x = (x±, ...,x n ), otherwise set 
k = k + 1 and repeat. 



Of course, in general the vector x € K n does not 
belong to K. Therefore a final step consists of computing 
a local minimum x G K, by using some local minimization 
algorithm starting with the (unfeasible) initial point x. Also 
note that when K is not convex, the determination of bounds 
x k and Xk for the interval Yfe may not be easy, and so 
one might be forced to use a subinterval Y' k C Yfe with 
conservative (but computable) bounds x' fe > x k and x' k < Xk- 

Remark 1: Theorem [2] assumes that for every y £ Yfe, 
the set K. y in (13[ is not empty, which is the case if K is 
connected. If K y = for y in some open subset of Yfe, then 
the semidefinite relaxation (fl4l > has no solution (pik = +oc), 
in which case one proceeds by dichotomy on the interval Yfe 
until pik < oo. 

C. A "joint+marginal" algorithm when K is convex 

In this section, we now assume that the feasible set K C 
l n of problem P is convex (and compact). The idea is to 
compute X\ as in ALGO 1 and then repeat the procedure 
but now for the (n — Invariable problem P(5i) which is 
problem P in which the variable X\ is fixed at the value x\. 
This alternative is guaranteed to work if K is convex (but 
not always if K is not convex). 

For every j > 2, denote by Xj £ M™ - ^ 1 the vector 
(xj, . . . , x n ), and by x,-_i G M J_1 the vector {x\, . . . , 
(and so Xi = x\). 

Let the interval Yi C 1 be the orthogonal projection of 
K into the x\ -coordinate axis. For every x\ G Yi, let the 
interval Y2(xi) cKbe the orthogonal projection of the set 
K n {x : x\ — x\] into the X2 -coordinate axis. Similarly, 
given x 2 G Yi x Y 2 (xi), let the interval Y 3 (x 2 ) cKbe the 
orthogonal projection of the set Kn{x : x± = X\; x% = a^} 
into the X3-coordinate axis, and etc. in the obvious way. 

For every k = 2, . . . ,n, and Xfe_i G Yi x Y2(xi) • ■ • x 
Yfe_i(xfe_ 2 ), let /fe(xfe) := /((x fe _i, x fc )), and g k (x k ) ■= 
3j((xfe_i,Xfe)), j = 1, . . . ,to. Similarly, let 

Kfe(x fe _i) := {xfe : (x fe ) > 0, j = 1, . . . , m}, 

= {xfc : (xfe_ 1)X fe) G K}, (18) 

and consider the problem: 

P(xfe_i) : min{/fe(x x ) : x K G Kj(x fc _i)}, (19) 

i.e. the original problem P where the variable X£ is fixed at 
the value xg, for every t = 1, . . . , k — 1. 

Write Yj(xfe_i) = [x k ,Xk], and let ipk be the probability 
measure uniformly distributed on Yfc(xfc_i). 

Let z be a sequence indexed in the monomial basis 
of M[xfc]. With index i, fixed, the parametric semidefinite 
relaxation ( Tl4l > with parameter Xk, associated with problem 
P(x fe _i), reads: 

p lk = inf L z (f k ) 

si Mi (z) , Mi- Vj (g k z)h0, j = 1 , . . . , to 
L z (4)=&, £=0,l,...,2i, 

(20) 



where (f3i) is defined in ( fTTT i. Its dual is the semidefinite 
program (with gfi = 1)): 

Pik = SU P / Pidtfk (21) 

Pi,(<?j) JY fc (x fc _i) 

m 

s.t. / fc - = cr + ^ o-j g| 

i=i 

Pi G M[£c fc ] 2 i, crj- G S[xfc], j = 0,...,m 
dego-js} < 2i, j = 0,...,m. 

The important difference between (fT4l i and (f20b is the s;ze 
of the corresponding semidefinite programs, since z in (O 
(resp. in (l20l i) is indexed in the canonical basis of M[x] (resp. 
K[xfc]). 

77ie (J+M)-algorithm for K convex 

Recall that the order i of the semidefinite relaxation is 
fxed. The (J+M)-algorithm consists of n steps. At step fc of 
the algorithm, the vector Xfe_i = (xi, . . . ,Xk-i) (already 
computed) is such that x% G Yi and xi G Y^(x£_i) for 
every I — 2,...,k — 1, and so the set Kfc(xfc_i) is a 
nonempty compact convex set. 

ALGO 2: (J+M)-algorithm: convex K, relaxation i 

Set k = 1; 

Step fc > 1: Input: For k = 1, x = 0, Yi(xo) = Yi; 
P(x ) = P, /i = / and <?j = 5 3 , j = 1, . . . , m. 
For k > 2, x fc _x G Yx x Y 2 (xi) ••■ X Y fc _i(x fe _ 2 ). 
Output: x fe = (x fc _i,£fc) with £ fc G Y A .(x fc _ 1 ). 
Consider the parametric semidefinite relaxations (l20T>-(l2Tb 
with parameter associated with problem P(xfe_i) in $1% . 

• From an optimal solution of (l2T1 i. extract the univariate 
polynomial x k ^ p,(x k ) := Y,f=o X *e x k- 

• Get a global minimizer x k of Pi on the interval 
Y fc (x fe _i) = [x fc ,x fe ], and set x fc := (x fc _i,5 fc ). 

If fc = n stop and ouput x G K, otherwise set k — k + 1 
and repeat. 

As K is convex, x G K and one may stop. A refinement 
is to now use x as the initial guess of a local minimization 
algorithm to obtain a local minimizer x G K of P. In view 
of Theorem 12 the larger the index i of the relaxations (f2Qb - 
(fSTJ, the better the values /(x) and /(x). 

Of course, ALGO 2 can also be used when K is not 
convex. However, it may happen that at some stage fc, the 
semidefinite relaxation (120b may be infeasible because J k (y) 
is infinite for some values of y G Yfe(xfc_i). This is because 
the feasible set K(xfc_i) in ( TT~8T > may be disconnected. 

IV. Computational experiments 

We report on preliminary computational experiments on 
some non convex NP-hard optimization problems. We have 
tested the algorithms on a set of difficult global optimization 
problems taken from Floudas et al. [1]. To solve the semidef- 
inite programs involved in ALGO 1 and in ALGO 2, we 
have used the GloptiPoly software [3] that implements the 
hierarchy of semidefinite relaxations defined in [4, (4.5)]. 
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TABLE I 
ALGO 2 FOR CONVEX SET K 



A. ALGO 2 for convex set K 

Those problems are taken from [1, §2]. The set K is a 
convex polytope and the function / is a nonconvex quadratic 
polynomial x h-> x'Qx+b'x for some real symmetric matrix 
Q and vector b. In Table I one displays the problem name, 
the number n of variables, the number m of constraints, the 
gobal optimum /*, the index i of the semidefinite relaxation 
in ALGO 2, the optimal value obtained using the output of 
ALGO 2 as initial guess in a local minimization algorithm 
of the MATLAB toolbox, and the associated relative error. 
As recommended in Gloptipoly [3] for numerical stability 
and precision, the problem data have been rescaled to obtain 
a polytope contained in the box [—1, 1]™. As one may see, 
and excepted for problem 2.8C5, the relative error is very 
small. For the last problem the relative error (about 11%) 
is relatively high despite enforcing some extra upper and 
lower bounds Xi < xi < Ti, after reading the optimal 
solution. However, using x G K as initial guess of the 
local minimization algorithm in MATLAB, one still finds 
the optimal value /*. 

B. ALGO 1 for non convex set K 

Again in Table II below, n (resp. m) stands for the number 
of variables (resp. constraints), and the value displayed in 
the "ALGO 1" column is obtained in running a local 
minimization algorithm of the MATLAB toolbox with the 
output x of ALGO 1 as initial guess. 

In Problems 3.2, 3.3 and 3.4 from Floudas et al. [1, §3], 
one has 2n linear bound constraints and additional linear and 
non convex quadratic constraints. As one may see, the results 
displayed in Table II are very good. 

For the Haverly Pooling problem 5.2.2 in [1, §5] with three 
different data sets, one has n = 9 and m = 24 constraints, 
among which 3 nonconvex bilinear constraints and 18 linear 
bound constraints < Xi < 500, i = 1, . . . , 9. In the first 
run of ALGO 1 we obtained bad results because the bounds 
are very loose and in the hierarchy of lower bounds (/£) 
in [4] that converge to /*, if on the one hand /| = /*, on 
the other hand the lower bound /* < /* is loose. In such a 
case, and in view of the rationale behind the "joint+marginal" 
approach, it is illusory to obtain good results with ALGO 
1 or ALGO 2. Therefore, from the optimal solution x* in 
[1], and when < x* < 500, we have generated stronger 



Prob 


n 


TO 


/* 


% 


ALGO 1 


rel. error 


3.2 


8 


22 


7049 


1 


7049 


0% 


3.3 


5 


16 


-30665 


1 


-30665 


0% 


3.4 


6 


18 


-310 


1 


-298 


3.8% 


5.2.2 (1) 


9 


24 


400 


1 


400 


0% 


5.2.2 (2) 


9 


24 


600 




600 


0% 


5.2.3 (3) 


9 


24 


750 




750 


0% 


5.2.4 


9 


24 


750 




750 


0% 


7.2.2 


6 


17 


-0.3746 




-0.3746 


0% 


7.2.6 


3 


7 


-83.254 




-82.3775 


1% 



TABLE II 
ALGO 1 FOR NON CONVEX SET K 



bounds QAx* < Xi < 1.6x*. In this case, /* is much closer 
to /* and we obtain the global minimum /* with ALGO 
1 followed by the local minimization subroutine; see Table 
II. Importantly, in ALGO 1, and before running the local 
optimization subroutine, one ends up with a non feasible 
point x. Moreover, we had to sometimes use the dichotomy 
procedure of Remark Q] because if Y/- is large, one may have 
K y = for y in some open subintervals of Yfc. 

Problem 7.2.2 has 13 linear constraints and 4 nonlinear 
constraints with bilinear terms. To handle the non-polynomial 
function x^ 5 , one uses the lifting uf = Xi, Ui > 0, 
i = 5,6. Problem 7.2.6 has only 3 variables, 6 linear 
bound constraints, and one highly nonlinear constraint (and 
criterion). Here one uses the lifting v?X2 = 1, u > 0, to 
handle the term x^ 1 . Again one obtains the optimal value /* 
with ALGO 1 followed by a local optimization subroutine. 

C. ALGO 2 for MAXCUT 

Finally we have tested ALGO 2 on the famous NP-hard 
discrete optimization problem MAXCUT, which consists of 
minimizing a quadratic form x i-)- x'Qx on { — 1, 1}™, for 
some real symmetric matrix Q £ R raxn . In this case, Y/. = 
{ — 1,1} and the marginal constraint L z (x e k ) — in (l20b 
need only be imposed for 1 = 1, because of the constraints 
x\ = 1 for every k = 1, . . . , n. Accordingly, in an optimal 
solution of the dual (1211 , pi e R[xfc] is an affine polynomial 
Xk Pi(xk) = Ao+Ai^fe for some scalars Ao, Ai. Therefore 
after solving (EH one decides Xk = —1 if p%(— 1) < 
(i.e. if Ai > 0) and Xk = 1 otherwise. 

Recall that in ALGO 2 one first compute x\, then with x\ 
fixed at the value x\, one computes x%, etc. until one finally 
computes x n , and get x. In what we call the "max-gap" 
variant of ALGO 2, one first solves n programs (fl4]i-(fT5Tl 
with parameter x\ to obtain an optimal solution Pi(x\) = 
Aq + A}a;i of the dual (TT3T >. then with x 2 to obtain (Aq, Af), 



n 


20 


30 


40 




10.3% 


12.3% 


12.5% 



TABLE III 
Relative error for MAXCUT 



etc. finally with x n to obtain (Aq , A"). One then select k such 
that |A^| = m&xg \X[\, and compute Xk accordingly. This is 
because the larger |Ai|, (i.e. the larger \pi(— 1) —pi(l)\), the 
more likely the choice —1 or 1 is correct. After Xk is fixed 
at the value ik, one repeats the procedure for the (n — 1)- 
problem P(ifc), etc. 

We have tested the "max-gap" variant for MAXCUT 
problems on random graphs with n = 20, 30 and 40 nodes. 
For each value of n, we have solved 50 randomly generated 
problems and 100 for n — 40. The probability ipk on Y^ = 
{ — 1,1} is uniform (i.e., /3i = in (f20b). Let /-j" denote the 
optimal value of the Shor's relaxation with famous Goemans 
and Williamson's 0.878 performance guarantee. Let p denote 
the cost of the solution x G { — 1,1}" generated by the 
ALGO 2. In Table III we have reported the average relative 
error (p — /* )/|/i |, which as one may see, is comparable 
with the Goemans and Williamson (GW) ratio. 

V. Conclusion 

First preliminary results are promising, even with small 
relaxation order i. When the feasible set is non convex, it 
may become difficult to obtain a feasible solution and an 
interesting issue for further investigation is how to proceed 
when K y = for y in some open subinterval of Y^ 
(proceeding by dichotomy on Y^ is one possiblity). 
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